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Abstract 

We present the software library libCreme which we have previously used to successfully calculate convex-roof entanglement 

measures of mixed quantum states appearing in realistic physical systems. Evaluating the amount of entanglement in such states 

is in general a non-trivial task requiring to solve a highly non-linear complex optimization problem. The algorithms provided here 

are able to achieve to do this for a large and important class of entanglement measures. The library is mostly written in the Matlab 

^^ programming language, but is fully compatible to the free and open-source Octave platform. Some inefficient subroutines are 

f~^ written in C/C-h-h for better performance. This manuscript discusses the most important theoretical concepts and workings of the 

p<| algorithms, focussing on the actual implementation and usage within the library. Detailed examples in the end should make it easy 

,__i for the user to apply libCreme to specific problems. 
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1. Introduction 

The role of non-local quantum correlations, more familiarly 
known as entanglement, in modern quantum theory cannot be 
overstated lU. On the one hand, entanglement lies at the heart 
of quantum information theory |T|, where it is a crucial ingredi- 
ent to computation and communication schemes. On the other 
hand, it is intricately related to phenomena such as decoher- 
ence |3l and quantum phase transitions in many-body systems 



|4|. One has come to realize that entanglement is also a re- 
source that can for instance be purified, shared, and possibly 
irreversibly lost, and should therefore not only be detectable, 
but also quantifiable (51. One way of doing so is by virtue of 
entanglement measures ||6l. These are mathematical functions 
mapping quantum states to the set of real numbers. While there 
is no unique or strict definition for the notion of an entangle- 
ment measure, there are a set of properties which are commonly 
regarded useful, e.g., that the measure is zero only for separa- 
ble states and is invariant under local unitary transformations. 
Another important property which we will assume throughout 
this work is monotonicity: An entanglement measure must not 
increase (on average) under any protocol involving only local 
unitary transformations and classical communication. In the 
following, we will use the terms 'entanglement measure' and 
'entanglement monotone' interchangeably. 

Rather understandably, it is difficult to capture all proper- 
ties of even a pure entangled state with just a single real num- 
ber, especially in the setting of higher-dimensional and multi- 
partite systems. It is thus no surprise that there is quite a number 
of proposed entanglement monotones of various levels of com- 
plexity, generality, and the ability to capture different aspects of 
entangled states more or less successfully than others. As indi- 
cated previously, most of these entanglement monotones share 
the fact that they are conveniently defined only for pure states, 
namely as a function of the amplitudes of the state expressed in 
a certain standard basis. 

The situation becomes more involved in the case of mixed 
states, where classical and quantum correlations need to be dis- 
tinguished from one another. Given a density matrix p, it is 
not sufficient to simply calculate the average entanglement of a 
given decomposition, because this decomposition is not unique. 
Since there are in general infinitely many ways to write a den- 
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sity matrix as a sum of projectors onto pure states, only the 
infimum of entanglement over all these possible decomposi- 
tions can make a reliable statement about the quantum corre- 
lations in p, because there might be a decomposition of p in 
which all pure states are separable and the total entanglement 
hence would vanish. Taking this infimum of an averaged pure- 
state entanglement monotone over all decompositions of p is 
called 'convex-roof construction' or 'convex-roof extension' of 
that monotone |7|. Notably, the thereby obtained measure for 
mixed states is again an entanglement monotone. Calculating 
the convex-roof for a generic quantum state is considered ex- 
tremely difficult f6l. In fact, even deciding whether or not a 
merely bipartite mixed state is separable is a hard problem itself 
which has no known general solution in Hilbert space dimen- 
sions larger than six fSl. 

In this work, we present the computer programs we have 
written and successfully applied previously to calculate such 
convex-roof entanglement measures of multi-partite mixed states 
IHj i9j . While we have already described the theory behind our 
algorithms to some extent in an earlier publication |l9], we com- 
plete this work by publishing here the full source code in the 
form a user-friendly high-level hbrary called libCreme. The 
package is to a large part written in the Matlab programming 
language fTOl, but great care has been taken to make the library 
fully compatible with GNU Octave | 1 1 |, a free and open-source 
Matlab clone. For the sake of simplicity, we will refer to their 
common language as M-script. Additionally, functions which 
have been identified as crucial bottlenecks in terms of execu- 
tion speed are provided in the form of fast C extensions and 
have been adapted to be easily callable from Matlab and Oc- 
tave through their native C and C-I-+ interfaces, respectively. 

While the library already comes with the abiUty to evaluate 
a choice of popular entanglement monotones, it is easily extend 
to calculate user-specified functions. All that needs to be imple- 
mented is the entanglement measure itself and its gradient with 
respect to the real and imaginary parts of the quantum state vec- 
tor |12 |. The library is written in a self-contained and consistent 
way, making it extremely easy to use in practice and to exper- 
iment with different settings, measures, and optimization algo- 
rithms. Furthermore, we provide convenience functions hiding 
most of the steps required to arrive at function handles ready 
to be optimized. This essentially simplifies the calculation of a 
convex-roof entanglement measure to a one-line task. Table [T] 
lists each function provided in the library together with a short 
description of its meaning. 

We would briefly like to mention two other numerical li- 
braries dealing with quantum computing and entanglement. One 
is the freely available 'quantum information package' by T. Cu- 
bitt |fT3]| written in M-script as well. The other one is the CPC 
library Feynman (catalog identifier ADWE_v5_0) written by T. 
Radtke and S. Fritzsche |14| for Maple. Quantum states ob- 
tained from calculations and simulations within these libraries 
can conveniently be anaylized further using libCreme's ability 
to calculate more complex entanglement measures. 

The paper is organized as follows: In Sec.l2] we briefly list 
and discuss the default entanglement measures coming along 
with libCreme. Sec. [3] reviews the theory of convex-roof en- 



tanglement measures and how to address their calculation on 
a computer Sec. |4] describes the two central algorithms in 
libCreme to solve the optimization problem related to the eval- 
uation of such measures. In Sec. B] we discuss two complete 
examples demonstrating the usage of the library, and Sec. |6] 
concludes the work. Note that the focus of this manuscript lies 
mainly on the functionality of the library: We have tried to pro- 
vide short code examples throughout the work for all important 
functions and concepts in a user-friendly bottom-up way. These 
snippets are all valid M-script (including the line breaks ' . . . ' 
which we sometimes use due to spacial restrictions) and can 
be copied directly into Matlab or Octave. Finally, we would 
like to mention that all functions in libCreme are documented, 
and more information about them can be inquired by calling 
'help function Jiame' . 

2. Entanglement measures included in the library 

We start the description of our library by listing the en- 
tanglement measures currently implemented. Note that pure 
quantum states, such as the arguments of functions calculat- 
ing entanglement monotones, are always expected to be rep- 
resented as column vectors in the standard computational ba- 
sis. In practice, this means that the n orthonormal basis states 
\ifri) of each n-dimensional subsystem (where n may be dif- 
ferent for different subsystems) are always chosen as \tfri} = 
(1,0,...,0)M<A2) = (0,1,0,. ..,0)^...,|<A„) = (0,0,...,!)^. 
Multipartite states are then assumed to be represented consis- 
tently with the implementation of the M-script command kron, 
i.e., the Kronecker product of two arbitrary input matrices. 

Since the optimization algorithms used in libCreme are 
gradient-based, the gradients of these measures with respect 
to the real and imaginary parts of the input state vector are 
also provided. The convention is that gradients (i) are named 
identical to the original functions but with the prefix 'grad_' 
added, (ii) require the same arguments as their function counter- 
parts, and (iii) return derivatives with respect to real and imag- 
inary parts of a variable in the form [V/(ji:)],- = df/dRexi + 
i df/d Im Xj, where i is the imaginary unit. Analytical expres- 
sions for all gradients of the measures discussed in this section 



can be found in Appendix A 



2.1. Entropy of entanglement 

The entropy of entanglement ifTSl is an entanglement mono- 
tone for bipartite quantum systems of arbitrary dimensions. It 
is defined as the von Neumann-entropy of the reduced density 
matrix, i.e.. 



£(|^))--Tr[(TrBp)log2(TrBp)]. 



(1) 



where Tr^p denotes the partial trace of p = |iA)('AI over the sec- 
ond subsystem (note that one could equally use the trace over 
the first subsystem). This measure is implemented in 
entropyOf EntcLnglement and requires as a first argument the 
state vector to be evaluated, and as the second a two-dimensional 
row vector specifying the dimensions of the two subsystems. 
Note that entropyOf EntcLnglement makes use of pTrace, a 



Entanglement measures 

convexSum 

grad_convexSum 

eof2x2 

entropyOf Entanglement 

grad_entropyDf Entanglement 

meyer_wallach 

grad_meyer_wallach 

tangle 

grad_tangle 



Convex sum parameterized by a Stiefel matrix 

Gradient of convex sum 

Entanglement of formation for 2 qubits (analytically exact result) 

Entropy of entanglement 

Gradient of entropy of entanglement 

Meyer- Wallach measure 

Gradient of Meyer- Wallach measure 

Tangle 

Gradient of Tangle 



Optimization routines 

cgjnin 

bf gsjnin 

minimize ld_exp 

minimize ld_lin 

get _termination_cr iter ia 



Conjugate-gradient method 

BEGS quasi-Newton method 

Minimization along a geodesic on the Stiefel manifold 

Minimization along a line in Euclidean space 

Helper function to check and handle termination criteria for the optimization 

algorithms 



Utilities 

raiidDensityMatrix 

raiidState 

randUnitaryMatrix 

decomposeUnitary 

dimSt 

densityEig 

psDecomposition 
createConvexFunctions 
createEHFunctions 
grad_eh_adapt 

buildUnitary 
grad_buildUnitary 
pTrace 
completeGrcunSchmidt 



Random density matrix 
Random pure quantum state 
Random Stiefel matrix 
Get angles parameterizing a Stiefel matrix 
Dimension of Stiefel manifold 

Get eigendecomposition of a density matrix in the form required by many func- 
tions within the library 

Get pure-state decomposition parameterized by a Stiefel matrix 
Create convex-sum function handles for use with cgjnin 
Create convex-sum function handles for use with bf gs_min 
Adapter function to calculate the gradient of a convex sum parameterized by 
Stiefel matrix angles 

Build a complex Stiefel matrix from a parameterization vector 
Gradient of the above function 

Partial trace over any subsystems of arbitrary (finite) dimensions 
Helper function for numerical stability used within cg_min 



Examples 

example_eof Isotropic 

eof Isotropic 

example_tangleGHZW 
tangleGHZW 



Main script to run the example from Sec. 5. 1 

Entanglement of formation of an 'isotropic density matrix (analytically exact 

result) 



Main script to run the example from Sec. 5.2 



Tangle of GHZ/W mixture (analytically exact result) 



Table 1: List of all functions within libCreme. Additional information about the usage of each function can be obtained by calling 'help function_name' from 
within Matlab or Octave. 



C/C++ implementation for the fast calculation of partial traces 
over an arbitrary set and number of subsystems of arbitrary di- 
mensions. Usage: 



Create states p_01 and p_10 in the total 
Hilbert space of two qubits by applying the 
Kronecker product to the single-qubit basis 
states [1; 0] and [0; 1] . 



Note that in M-script, [al ; a2; . 
denotes a column vector, whereas 
[al , a2, ... an] is a row vector. 
,01 = kron([l; 0], [0; 1]); 
.10 = kron([0; 1], [1; 0]); 



. an] 



7o Define a random superposition of the 

°/o above states. 

7. 

% randO yields a random number chosen 

% uniformly from the interval (0, 1). 

rl = 2*pi*rand(); r2 = 2*pi*rand() ; 

psi = sin(rl)*p_01 + exp(li*r2)*cos(rl)*p_10; 

7o Dimensions of subsystems 
sys = [2, 2]; 

7o Calulate measure and gradient 

e = entropyOf Entanglement (psi, sys) 

g = grad_entropyOf Entanglement (psi, sys) 

Note that the entropy of entanglement is of particular impor- 
tance, because its convex-roof extension is the well-known and 
widely used 'entanglement of formation' |16|. In the special 
case of a bipartite system composed of two-dimensional sub- 
systems (qubits), there exists an operational solution for the en- 
tanglement of formation 1,17,1 . which we have implemented in 
eof2x2. 

2.2. Three-tangle 

The three-tangle IfTSl is defined specifically for a system of 
three two-dimensional subsystems. It reads 



where 



d\ — 
d2 = 

di - 



Tm)^4\di-2d2+4di\, (2) 



+lA4lA5lA6'A3 + >p4^5^1<p2 + >Pb>Pi^l^2, (4) 

4ll4l-lil/(,ll/A + iA8l/'2'A3lA5, (5) 



and 4'i,i - 1 • . . , 8, are the complex amplitudes of the vector |i/') 
in the standard computational basis. The tangle is implemented 
in the function teingle, taking an 8-dimensional vector psi as 
its only argument. Here follows a short example: 

% Define an 8-dimensional reuidom state 
psi = randState(8) ; 



7. Calculate tangle and its gradient 

t = tangle (psi) 

g = grad_tEuigle(psi) 

In the above code, we have introduced the function randState, 
which returns a random pure quantum state of arbitrary speci- 
fied dimension uniformly distributed according to the Haar mea- 
sure of the unitary group |fT9l . 

2.3. Meyer-Wallach measure 

Finally, the measure of Meyer and Wallach fTO] for an arbi- 
trary number A^ of qubits is an entanglement monotone that can 
be written in the compact form 11211 



rm) = 2 



l-^ZTr(p,^) 

k=l 



(6) 



where p^ is the density matrix obtained by tracing out all but 
the ^th subsystem out of \t//}{i//\. The implementation is given 
in meyer_wallach. This function also makes use of pTrace. 
The usage is analogous to the example given for the three-tangle 
above. 

3. Theoretical background 

In this section, we review how to arrive at an optimization 
problem (whose solution is the desired value of the convex-roof 
entanglement measure) in a form that can be dealt with on a 
computer. Let m be an entanglement monotone for pure states 
from a Hilbert space '7/ of finite dimension d. Let p be a density 
matrix acting on that space. Our goal is to numerically evaluate 
the convex roof M(p) of m, given by 



M(n) = inf 

|pi,|^,»eS(p) 



Y^PimQiffi)), 



(7) 



where 



S)(P) = [{piA'Pi}]li,k> rankp | {|,A/)ti c 'H, 

k k 

{i/,i\i/,i} = I, Pi > 0, ^ Pi = 1, p = ^ P/l-A/X-A/l} (8) 
(=1 1=1 

is the set of all pure-state decompositions of p. With respect 
to numerical optimization, a convenient parameterization of all 
subsets of ©(p) with a constant number of terms k (sometimes 
referred to as the 'cardinality') is due to the Schrodinger-HJW 
theorem Il22ll23l . The latter states that (i), every decomposi- 
tion of a density matrix p with rankp = r into a convex sum 
of k projectors onto pure states can be expressed in terms of 
a complex k x r matrix U obeying U^U - Irxr and that (ii), 
conversely, from every such matrix one can obtain a pure-state 
decomposition of p. The set St{k,r) = {t/ e C*^nt/"'"t/ = I,x.} 
with A: > r is also known as the Stiefel manifold. Part (i) and 
(ii) together ensure that optimizing over S t(k, r) is equivalent to 
optimizing over the full subset of 2)(p) with fixed cardinality k. 
Part (ii) also provides an explicit construction of the pure-state 



decomposition related to an arbitrary given matrix U € St{k, r): 
Let Ai, \xi), i = 1, . . . , r = rankp be the non-zero eigenvalues 
and corresponding normalized eigenvectors of p, i.e., 



P = ^ ^i\Xi){Xil 



(9) 



and ixilxj) - Sij. Note that we have /I, > due to the positive 
semi-definiteness of p. Define the auxiliary states 

\ijfi)^Y^Uij4Aj\xj), i^l,...,k, (10) 



and set 

Pi = {>Pi\>Pi), 

\^d = (1/Vp^)I'A/>- 

One can easily verify that we have 



P^^Pi\^i){^i\- 



(11) 
(12) 



(13) 



(=1 



In libCreme, the function densityEig calculates only the 
non-zero eigenvalues and corresponding eigenvectors. The eigen- 
values are guaranteed to be sorted in decreasing order, which is 
particularly convenient if one wishes to discard some parts of 
the density matrix occurring with low probability, such as, e.g., 
high-energy sectors in density matrices p ~ exp{-/3H) originat- 
ing from some Hamiltonian H. The function psDecomposition 
returns the pure-state decomposition from Eqs. ([9] [T0l[TT][T2| . 
As an example, let rho store a d x d density matrix of rank r, 
and let U be a matrix from S t(k, r), with arbitrary k> r. Then 

°/o Note that in M-script, functions CEin return 

7o multiple values of arbitrary dimensions. The 

7o syntax to assign several return values to 

7o local variables is 

I [A, B, ...] = foo(...); 

[chi, lambda] = densityEig(rho) ; 

[psi, p] = psDecompositionCU, chi, lambda); 

first yields the eigenvectors of rho in the columns of the dx r 
matrix chi with the corresponding r eigenvalues in the vector 
lambda. On the second line then, the pure-state decomposition 
of rho (given in terms of the parameters chi and lambda) cor- 
responding to the parameterization U is calculated, with the k 
state vectors \\jjj) stored in the columns of the dxk matrix psi 
and the k corresponding probabilities pi in the vector p. 

By virtue of the Schrodinger-HJW theorem, we can restate 
the optimization problem Eq. (|7]i as 



Mip) - min inf hiU), 

k>r UeSt(k,r) 



h(U) = YjPi(U)m(\ifri(U))), 



(14) 
(15) 



practice, it has turned out to be possible to drop the minimiza- 
tion over k completely and set ^ to a constant but large enough 
value instead. Note that this actually includes all cardinalities 
k' with r < k' < k in the search because up to ^ - r of the 
Pi are allowed to go to zero without breaking the optimization 
constraint. In libCreme, the function convexSum calculates 
the value of the expression h{U) in Eq. ( [T5] l, which is, in fact, 
the objective function of the optimization. convexSum takes 
as its first argument a parameterization matrix from the Stiefel 
manifold, as its second a function handle I.24J to the entangle- 
ment monotone to be extended, and as its third and fourth ar- 
guments the eigendecomposition of the density matrix obtained 
by densityEig. Here is a full example: 

7„ Reindom 8-by-8 density matrix 
rho = randDensityMatrixCS) ; 

7o Calculate eigendecomposition of rho for 

7o later use 

[chi, Icmibda] = densityEig(rho) ; 

7» Rcindom matrix from St (12, 8) 
U = randUnitaryMatrix(12, 8); 

7. Evaluate convex sum Eq. (15) with the tcingle 
h = convexSum (U, Stangle, chi, lambda) 

Note that we have introduced the functions randDensityMatrix 
and randUnitaryMatrix to create random density matrices 
and random matrices from the Stiefel manifold, respectively. It 
is important to understand that convexSum is the key function 
in the whole library in the sense that it is always this function 
(or more specifically, an anonymous function handle |24| to it, 
see Sec. |4|, which is ultimately optimized. 

As mentioned earlier, the optimization algorithms in 
libCreme require the knowledge of the gradient of the objec- 
tive function, or more precisely, the derivatives of h(U) with 
respect to the real and imaginary parts of the matrix elements 
of U. These expressions and their derivation can be found 
Within the library, this gradient of h is im- 



Appendix B 



where the dependence on p enters implicitly as the dependence 
of the Pi and |i^,) on the eigenvalues and eigenvectors of p. In 



m 

plemented in grad_convexSum. It requires 5 arguments: The 
matrix U, the entanglement monotone to be extended, the gra- 
dient of the latter, and the eigendecomposition of p (eigenvec- 
tors and -values). The following code illustrates its application 
in practice, using the variables chi, lambda, and U from the 
previous example: 

7o Evaluate gradient of the convex sum Eq. (15) , 
7o given in Eqs . (B4 - B7) with the tcoigle 
gh = grad_convexSum(U, Stangle, ... 

@grad_tangle, chi, lambda) 

4. Optimization algorithms 

We describe in this section two conceptually diff'erent op- 
timization algorithms which are both provided in libCreme. 



One is a conjugate gradient method based on the concepts in- 
troduced in Refs. ||251 l26l l27l . It exploits the differential- 
geometric structure of the nonlinear search space emerging from 
the optimization constraint f/'f/ = I. The other algorithm is 
a standard Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi- 
Newton method employing a transformation of the constrained 
search space to an unconstrained one. Both algorithms have 
been discussed in greater detail in a previous work |9|, where 
the expressions for the gradients and parameterization of the 
Stiefel manifold given below have been derived. The interested 
reader is referred to that earlier work. Here, we just state the 
final results for the sake of completeness and focus particularly 
on the implementation and usage within libCreme. 

4.1. Conjugate-Gradient Method 

This algorithm exploits the geometric structure of the uni- 
tary group U{k) - S t{k, k) and therefore generally over-para- 
meterizes the true search space S t(k, r), r < k. This is however 
not a problem in practice, since we can simply discard the last 
k- r columns of U when calculating the decomposition of the 
density matrix based on U Il28ll . The full algorithm for an input 
initial guess t/o is given as follows: 

1. Initialization: Set ; < — 0. Calculate the gradient Go = 
G(Uo) according to the formula 

G(U)^^-(A(U)-A(Uf)+^(S(U) + S(Uf), (16) 



where the matrices A{U) and S(U) are given by 

A( f/) = Re f/^ ■ VRe uh(U) + Im f/^ -Vi^uKU), (1 7) 
S(U) = Re [/^ ■ Vi^uh(U) - Im t/^ ■ VRef//z(t/), (18) 



and the gradient of h{U) can be found in Appendix B 
Finally, set Xq < 



-Go. 
2. Perform the one-dimensional minimization 

ti+i < — argmin/i(C/,exp(fX,)), 
set 



(19) 



f/,+1 ^f/,exp(f,+iX,), (20) 

and compute the new gradient G,+i < — G(Ui+i) accord- 



ing to Eqs. ( [T6|[T7| [T8i. 
3. Define 



exp(f,+iX,/2)G,exp(-f/+iX,/2). 



(21) 



T is the gradient G, parallel-transported to the new point 
4. Calculate the modified Polak-Ribiere parameter 



r 



(Gi+i - T,Gi+i) 
{Gi,Gi) ' 



(22) 



where <X, 7) = TrXr'. 
5. Set the new search direction to 



X: 



i+\ 



-Gi^i + yXj. 



(23) 



6. Set/ < — / + 1. 

7. Repeat from step 2 until convergence. 

This algorithm is implemented in the function cgjmin. The 
minimization in step [2] is performed by the derivative-based 
method minimize ld_exp. cg_min requires the function to be 
minimized, its gradient, an initial point, and optionally a struct 
with user-specified termination criteria discussed below. At this 
point, we would like to work through a full example demon- 
strating the use of cg_min to calculate the convex-roof extended 
three-tangle of a mixed state. 

7, Random 8-by-8 density matrix 
rho = rsindDensityMatrixCS) ; 

7, Calculate eigendecomposition of rho for 

7. later use 

[chi, Icimbda] = densityEig(rho) ; 

7o Define anonymous function handles [24] to 
7o the objective function and its gradient 
f_opt = @(x) convexSumCx, Otangle, ... 

chi, lambda) ; 
g_opt = Q(x) grad_convexSum(x, Otangle, ... 

@grad_tangle, chi, Isimbda) ; 

7o Choose a rsindom starting point, here for 
7o a decomposition with cardinality 12 
UO = randUnitaryMatrix(12, 12); 

7o Perform the optimization 

[t , Ut , info] = cg_min(f _opt , g_opt, UO) ; 

A few comments about the above code are in order First, 
note that because cgjmin requires the objective and its gra- 
dient in the form of one-parameter functions, we need to de- 
fine the anonymous function handles f _opt and g-opt before 
continuing. In this way, f _opt is a new function that evalu- 
ates convexSum at a variable unitary input matrix while keep- 
ing the constant arguments Otaingle, chi, and lambda fixed. 
A similar description holds for g_opt. Second, note that the 
initial search point UO is a unitary and therefore square ma- 
trix, although matrices from 5'r(12, 8) would be sufficient in 
the example above to parameterize pure-state decompositions. 
The reason is, as mentioned above, that cgjnin is operating 
on the unitary group instead of the Stiefel manifold. However, 
this is hidden from the user by the fact that both convexSum 
and grad_convexSum can accept larger input than required, 
and automatically discard any dispensable columns. Third, we 
would like to draw the reader's attention to the output values 
of cgjnin. In the above example, t stores the convex-roof of 
the entanglement monotone to be evaluated (in this case the 
three-tangle) and Ut (or more precisely, the first r columns of 
it) represent the pure-state decomposition of rho arriving at this 
value. The variable info is a struct that carries useful additional 
information, namely the criterion that terminated the iteration 
(info . status), the function values along the iteration, exclud- 
ing intermediate values during line searches (info . f vals), and 



finally, the traversed points in the search space corresponding to 
fvals (inf o.xvals). 

An optional fourth argument containing settings for the ter- 
mination of the algorithm can be passed to cgjnin. The fol- 
lowing code illustrates the possible struct variables (the values 
in this example are also the defaults for any variables not set). 

7o Create a struct 
opts = struct ; 

7, Maximum number of iterations 
opts.MaxIter = 1000; 

°L Tol. on consecutive function values 
opts.TolFun = le-12; 

7, Tolerance on norm of difference between 
7o two consecutive gradients 
opts.TolG = le-10; 

7o Toleraince on norm of difference between 
7o two consecutive iteration points 
opts.TolX = le-10; 

The iteration is stopped if either the maximum number of iter- 
ations is reached, or one of the checked values is lower than its 
respective tolerance. Finally, we would like to mention that for 
the convenience of the user, there is a function called 
createConvexFunctions which performs all the necessary 
steps before the actual optimization in one line: 

7. Again using the tangle as an example 

[f_opt, g_opt] = createConvexFunctions (rho , ... 

Otcingle, @grad_tangle) ; 
[t , Ut , info] = cg_min(f _opt , g_opt, UO) ; 

As with any other numerical optimization procedure, it is advis- 
able to repeat the computations with different (random) initial 
conditions in order to reach a better approximation to the global 
minimum. 

4.2. BFGS Quasi-Newton Method 

The second algorithm is a classical BFGS quasi-Newton 
method |29| that makes use of a transformation which is able to 
unconstrain the optimization problem Eq. ^\A\ from the Stiefel 
manifold to ordinary Euclidean space. This transformation is 
conceptually identical to the example where one has an opti- 
mization problem with the constraint x^ +y^ - 1 and then sets 
X - sm0,y - cos 6 and performs the optimization over the new 
variable 6. Again, we only state in the following the main re- 
sults required to implement the algorithm and refer the reader 
interested in a thorough derivation to Ref. |9|. 

The number of independent real parameters required to pa- 
rameterize S t{k, r) is equal to its dimension which is given by 
d\vi\St{k,r) - 2kr - r^. Given a tuple of 'angles' X,, i = 
1 , . . . , dim S t(k, r), we relabel them in the following (arbitrary 
but fixed) way: ?? = [Xj] for / = 1, . . . , r[;t - (r -H l)/2], tp = {X,} 
for / = r[k -{r+ l)/2] H- 1, . . . , r{2k - r - 1), and x = Kl 



for/ = r(2k-r- 1)+ 1,. 
U e S t(k, r) according to 

U{X)=U{§,ip,X)^ 



,A\TiiSt{k,r). Then, we calculate 



*■-/ 



nn^*-/^--'^-.> 



i=l i=l 



R 



(24) 



where c,y = (/ - \){k - i/2) + j,Risakxr matrix with the 
only non-zero elements being /?„ - exp(i;^f,) for / = 1 , . . . , r, 
and the 'inverse Givens matrices' G^' are defined in terms of 
their matrix elements as 



[G-\§,<p)]ij = 



e "^ cos -&, if i - j = s, 

-e'^f sin ■&, if / = i, 7 = i -H 1 , 

e'^sint?, if / = s + l,j^ s, (25) 

e'^ cos )?, if i = s + l,j - s +\, 



'tj^ 



otherwise. 



In libCreme, calculating a Stiefel matrix from a vector of 
angles X by Eq. ^A\ is implemented in buildUnitary as a fast 
C/C-I-+ extension and is demonstrated in the following: 

7o Dimensions of the Stiefel manifold 
k = 10; r = 7; 

7o Random vector of angles with proper size 

7. (Uses dimSt for the dimension of the 

7o Stiefel manifold) . See footnote [30] . 

7. 

7. randnCm, n) yields sin m-by-n matrix of 

7o normally distributed random numbers. 

X = 2*pi*randn(l, dimStCk, r)); 

°L Finally, the Stiefel matrix 
U = buildUnitary (X, k, r) ; 

The derivatives of U(X) with respect to the angles X,- are imple- 
mented in the function grad_buildUnitary. The inverse op- 
eration of buildUnitary, namely obtaining the parameteriz- 
ing angles for a given matrix U can be performed by 

decomposeUnitary as in 

X = decomposeUnitary (U) ; 

This function is implemented in regular in M-script, because it 
is called only infrequently and thus is not time critical. 

We have now all the tools to describe the full BFGS quasi- 
Newton algorithm starting from an initial vector of angles Xq. 

1. Set / ^ 0, Ho ^ I, Go = Vx/i(f/(X))|x=x„, and S^ = 
-Go. Hq is the initial guess for the approximate Hessian, 
h is the convex sum Eq. ([T5|, and U(Xo) is the transfor- 
mation Eq. p4l ). 

2. Perform the line minimization 



and set 



ti+i <— aigmm h(U(Xi + tS i)) 



Xi+i < — Xi + ti+\S i. 



(26) 



(27) 



3. Compute the new gradient 

4. Update the approximate Hessian as IISTI 



(28) 



"--"-(-^l 



66^ 



(29) 



6'^y 



where the column vectors 6 and y are defined as 6 
X,+ i - Xi and y = G,+i - G,. 
5. Set the new search direction to 



I/+1 



-///+iG,-+i 



(30) 



6. Set/ < — /+ 1. 

7. Repeat from step 2 until convergence. 

The line minimization in step 2 is performed by minimizeld_lin, 
a subroutine that is conceptually identical to the function 
minimize ld_exp used above in the conjugate-gradient method. 
The full algorithm is implemented in bfgsjnin and its input 
and output parameters are identical to the ones in cgjnin. 
Hence, the descriptions in the previous section can be adapted 
analogously to bf gsjnin. However, the target function (and 
its gradient) look slightly different in the current case and are 
somewhat more cumbersome in terms of function handles, be- 
cause the additional intermediate transformation Eq. ^A\ needs 
to be incorporated. The following is a fully working example 
that should help to clarify this issue. 

7o RcOidom 8-by-8 density matrix 
rho = rEindDensityMatrix(8) ; 
[chi, lambda] = densityEig(rho) ; 

7o Convex-sum function handles for the tangle 
f_cr = (§(x) convexSum(x, OtsLngle, ... 

chi , Icimbda) ; 
g_cr = @(x) grad_convexSum(x, Otsingle, ... 

(§grad_tangle , chi , lambda) ; 

7o Dimensions of the Stiefel manifold 
r = rank (rho) ; 
k = r + 4; 

7, Objective function and gradient 

f_opt = S(x) f _cr (buildUnitary (x, k, r)); 

g_opt = @(x) grad_eh_adapt (x, k, r, g_cr) ; 

7o Choose a random starting point 
XO = 2*pi*randn(l, dimStCk, r)); 

7o Perform the optimization 

[t , Xt , info] = bf gs_min(f _opt , g_opt, XO) ; 

Notice the use of the auxiliary function grad_eh_adapt which 
calculates the gradient VxHUiX)) given the derivatives of h{U) 



with respect to the real and imaginary matrix elements of U. 
For the convenience of the user, there is a function that hides 
all the above steps just like in the case of the conjugate-gradient 
algorithm. It is called createEHFunctions and is analogously 
called, as exemplified here: 

[f_opt, g_opt] = createEHFunctions (rho, ... 
k, r, Stangle, @grad_t£m.gle) ; 
[t , Xt , info] = bf gs_min(f _opt , g_opt, XO) ; 

Clearly, the same note as in the previous section regarding mul- 
tiple restarts holds. 

Finally, we would like to make a remark about a detail in 
our implementation of bf gsjnin. It has shown to be useful in 
practice to reset the angles modulo In every few iterations. This 
improves numerical stability and convergence in the vicinity of 
a minimum. It is also advisable to vary the interval size after 
which this is done as this can improve performance depending 
on the problem. If bf gs_min is to be employed for non-periodic 
functions, these lines of code must be removed. 

4.3. General Remarks 

Before we end this section and look at some more exam- 
ples, we would like to make a few comments. At this point, 
the reader might wonder why we provide our own implementa- 
tion of a line search and a BFGS quasi-Newton method, instead 
of resorting to the functions available in Matlab and Octave. 
There are several reasons for that. First of all, it makes the li- 
brary independent of the platform, since the standard routines in 
Matlab and Octave work differently and hence generally pro- 
duce unequal results. Furthermore, having a simple implemen- 
tation at hand allows the user to quickly make modifications 
and customize the code to specific needs. In Octave this is can 
only be achieved with quite an effort, whereas in Matlab it is 
generally not possible at all. Additional issues are availabil- 
ity and backward-compatibility. While Matlab's optimization 
framework is well established, it is only available through the 
purchase of the 'Matlab Optimization Toolbox'. On the other 
hand, there is a free octave package for non-linear optimization 
tasks ||32l . But since this is still under active development, its 
usage within libCreme might potentially become incompatible 
with future releases of the package. 

Next we would like to address the performance of the algo- 
rithms as a function of r = rankp and k > r. The dimension 
of the Stiefel manifold (and hence the problem size) is given 
by dim5'f(^, r) - 2kr - r^. Since we must have k > r, we 
can replace it by A; = r -H n, with « > 0, yielding dim S t{k, r) - 
r^+2nr. This shows that the computational cost grows quadrati- 
caUy with the rank of p, but only linearly with the (user-specified) 
cardinality. The algorithms in the library are thus most efficient 
for low-rank density matrices, whereas experimenting with dif- 
ferent cardinalities is not that costly. Actually, already choices 
for n as low as n ^ 4 have shown to produce very accurate 
results in practice (see also below). Since, on the other hand, 
the scaling with r is less favorable, it is advisable to examine 
whether the rank of r can be reduced. Particularly in density 
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Figure 1 : Comparison of numerical experiments with theory. The solid line 
demonstrates the convergence of the entanglement of formation of an isotropic 
state (example in Sec. |5.1( , whereas the dashed line does the same for the tangle 
of a GHZ/W mixture (example in Sec.|5.2|. 



matrices originating from physical Hamiltonians it is often jus- 
tified to neglect high-energy states. As mentioned earlier, re- 
ducing the rank of p can conveniently be achieved by truncating 
the output of the function densityEig. 

5. Examples 

In this section, we demonstrate the usage of libCreme by 
working through two complete examples. We calculate the en- 
tanglement of special states where analytical results are known 
in order to compare the numerical experiments with theory. 
Note that we provide initial points for the optimization in sep- 
arate files instead of generating them randomly, because (i) the 
random number generators in Matlab and Octave produce dif- 
ferent sequences of random numbers and (ii) not every initial 
point leads to the convergence to a global minimum in such 
high-dimensional spaces. 

5.7. Entanglement offormation of isotropic states using cg_min 

Isotropic states are defined as a convex mixture of a maxi- 
mally entangled state and the maximally mixed state in a system 
of two li-dimensional subsystems. The isotropic state with an 
amount of mixing specified by /, where < / < 1, is given by 

Ea 

1-/ 



Pf 



(f-\ 



{i-m{r\)+m^){^% 



(31) 



where \ip*) 



' Yfi=\ I")- An analytical solution for the entan- 



V5 



glement of formation as a function of / and d has been found 
||33]| and is implemented in eof Isotropic. Let us compare 
now the numerical results with theory. The full example can be 
found in the folder exsunples/eof Isotropic, along with all 
other related files. 

We first choose a dimension for the two subsystems, 

d = 5; 

then create the maximally entangled state [i/f^) in these systems 
and store it in psi: 



psi = 0; 



for i = l:d 



tmp = zeros(d, 1); 

tmp(i) = 1; 

psi = psi + kronCtmp, tmp); 



end 



psi = psi/sqrt(d); 

After choosing a value for the mixing parameter f , 

f = 0.3; 

we can construct the isotropic state specified by d and f as 

7, Note that in M-script, A' is the Hermitian 
7„ conjugate of A. 

rho = (1 - f)/(d-2 - !)*( eye(d-2) - ... 
(psi*psi') ) + f *(psi*psi ' ) ; 

and calculate its eigendecomposition with 

[chi, Icunbda] = densityEig(rho) ; 

In order to keep this example fully reproducible, we unfortu- 
nately have to load and overwrite the eigenvectors chi from a 
file at this point. The reason is that the density matrix is degen- 
erate, yielding different eigendecompositions for the degenerate 
subspace depending on whether one uses Matlab or Octave due 
to the different diagonalization routines employed by these plat- 
forms. Clearly, one arrives at comparable results in both cases, 
but the paths in optimization space are different. 

chi = load( 'example_eof Isotropic_chi . txt ' ) ; 

After setting an appropriate cardinality 

r = rank (rho) ; 
k = 2*r; 

and defining function handles for the entanglement measure and 
its gradient 

eoe = (§(x) entropyOfEntanglement (x, [d, d] ) ; 
grad_eoe = . . . 
@(x) grad_entropyOfEntEinglement (x, [d, d] ) ; 

we can create the function handles required in the optimization 

f_cr = @(x) convexSum(x, eoe, chi, Ismibda) ; 
g_cr = @(x) grad_convexSum(x, eoe, ... 

grad_eoe , chi , lambda) ; 

Finally, we choose a random initial value UO (here initialized 
from a file) 

UOr = load( 'example_eof Isotropic_UOr . txt ' ) ; 
UOi = loadC 'example_eof Isotropic_UOi . txt ' ) ; 
UO = UOr + li*UOi; 



and perform the optimization: 



[e_res, U_res, info] = cg_min(f _cr , g_cr, UO) ; 



6. Conclusions 



This yields a value of e_res ^ 0.129322085695260 after 80 
iterations. We can check the convergence and the accuracy of 
the result by plotting the difference between the function value 
in each iteration and the theoretical value: 

semilogy(abs(inf o. f vals - eof Isotropic (f , d))); 

This produces the solid line in Fig. 1, showing that the result is 
exact up to an absolute error of ^ 10"'^. 

5.2. Three-tangle of GH^W mixtures using bfgs_min 

In this example, we will calculate the three-tangle of a mix- 
ture of the two states 

|GHZ) = (1000) -H 111 1))/V2 (32) 

|W) = (1001) + 1010) + 1100))/ V3, (33) 

given by Il34l 

Pp = p|GHZ)(GHZ| + (1 - pWXWl (34) 

The example files are in examples/tangleGHZW. 
In the code, we define the states 

GHZ = [1; 0; 0; 0; 0; 0; 0; l]/sqrt(2); 
W= [0; 1; 1; 0; 1; 0; 0; 0]/sqrt(3); 

choose a particular value for p 

p = 0.7; 

and create the mixed state 

rho = p*GHZ*GHZ' + (1 - p)*W*W'; 

We then specify a value for the cardinality 

r = rank (rho) ; 
k = r + 4; 

and can create the objective function and gradient handles. Note 
that we use the auxiliary function createEHFunctions to do 
all the required work: 

[f_eh, g_eli] = createEHFunctions (rho, ... 

k, r, Otangle, (§grad_tangle) ; 

After choosing a random initial point (initialized from a file, as 
before) 

XO = load('example_tangleGHZW_XO.dat') ; 

we are ready to perform the optimization: 

[t_res, X_res, info] = ... 

bf gs_min(f _eh, g_eh, XO) ; 

The value one obtains in this way after 29 iterations is t_res 
~ 0.190667409058084. A comparison with the analytical value 
ll34l is exact within numerical precision. We can again plot the 
error between the function values and the exact result 

semilogy(abs(inf o. f vals - tcLngleGHZW(0.7) ) ) ; 

which yields the dashed line in Fig. 1. 



We have presented our library libCreme which serves to 
evaluate generic convex-roof entanglement measures. The li- 
brary contains all tools required to deal with this problem, in- 
cluding two optimization algorithms working on the space of 
density matrix decompositions. The first one is based on a 
conjugate gradient algorithm operating directly on the group 
of unitary parameterizations, while the second one is a stan- 
dard BFGS quasi-Newton method employed with a transforma- 
tion from the original search space to unconstrained Euclidean 
space. Both implementations accept generic function handles, 
making it easy to extend their application to user-defined entan- 
glement measures. All that needs to be done for this is the im- 
plementation of the respective pure-state entanglement mono- 
tone and the corresponding derivatives with respect to the real 
and imaginary parts of the input state vector 
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Appendix A. Derivatives of entanglement measures 

In the following, we provide the calculations for the deriva- 
tives of all entanglement measures included in libCreme. 

Appendix A.l. Entropy of entanglement 

Let I^F) be a state vector from a bipartite system with sub- 
system dimensions di and d2. Let us rewrite Eq. Q in the form 



where 



E{\^'))^S{1rBP), 



S(X) = -TrXlogX, 



(A.l) 



(A.2) 



and p - |*I')(^|. Let i/r be an arbitrary (complex) entry of the 
state vector |*P). Then, using the chain rule, we have 



di// 



i,j,k.l 



dS(X)) 



X: 



X=Tr BP 



d(Tr Bp)ijdpki 
dpki dtp 



(A.3) 



Note that the indices k and I in the above sum run over the full 
Hilbert space dimension d\d2, whereas i and j only run over the 
first subsystem with dimension di . We now evaluate each term 
in the sum separately. 
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For the gradient of S (X) we get 

VxS (X) = - Vx Tr [{X - I) log X + log X] 

n=l 



V ,7=1 

°^ /_i\n+l ^ 

+ 2^^Tr[(X-ir] 



+ ^(-ir'(x-i)"- 
«=i 



•logX' 



where we have made use of the formula VxTt{X") = «(X""')^ 
and the fact that the series expansion of the logarithm is valid 
because in our case X is always a density matrix, thus having 
real eigenvalues between and 1 . 

Next we evaluate the derivatives of the partial trace Tr^. We 
wiU write coordinate indices of vectors in the full Hilbert space 
as ^ = (ki - l)d2 + ki = {k\,k2\, where k\ e {\,...,d\] and 
k2 G {!,... ,d2]- Then 



di 
k=l 



(A.4) 



and thus 



d(Ti Bp)ij 

^Pmn 



— / , ^i,m] ^k,m2^ j,n\ ^k,t 
k=\ 



Finally, we have to consider the derivatives of the density 
matrix itself with respect to the entries of the state vectors. This 
is the part where we have to treat Re ip and Im i// as independent 
variables because p is not analytic in the entries of I^F). One 
quickly finds 



dpki 
d Re tp„ 

dpki 
d Im tp„ 






(A.5) 
(A.6) 



Putting all this together we find, after eliminating all Kro- 
necker 5-symbols, 



dReiff„ 



£ {[V5(p)]„„. ■ «A[,>,] + [V5(p)],„, ■ .A[,>,]} (A.7) 



!=1 



and analogously for the derivatives with respect to the imagi- 
nary parts. Exploiting the fact that V^CTr^p) - - logCTr^p)^ - 
I is Hermitian, we arrive at the final expressions 

^^ = 2 V Re {[VS(TrBp)U ' 'A[,>]1 > (A.8) 

^^ = 2£lm{[V5(TrBp)],„, ■ ^[,-,„,]) . (A.9) 

Appendix A.2. Three-tangle 

Defining d - d\ - Idi + 4-dy,, where d\,d2,dT, are given in 
Eqs. ([3l|4]|5]), it is easy to see that 



Re 



dT(\'¥)) _ 4 
<9Re(A„ ~ Rl 



dT(^)) ^_^^id±,^. 






d Im \j/n \d\ \ d(f/„ 



(A. 10) 
(A. 11) 



Note that the derivatives dd/dil/„ are well-defined because d is 
an analytic function of the elements (//„ of I*!*). 

Appendix A.3. Meyer-Wallach measure 

We start by calculating the derivative of yd*!*)) with respect 
to an arbitrary complex element ifr of I^P) until the point where 
non-analyticities appear: 



dip 



2 Af 2 „ 
k=\ j=\ 



2 " ^ (5 



N 
4 

'n 



k=\ j,l=\ 

N 2 



(A. 12) 



zz^H(4.i: 



'd(pk)ji dpap 



k=\ j.l=l 



afi=\ 



dpa/i dlff 



The derivatives of the paj^ (depending non-analytically on 
t//) with respect to the real and imaginary part of ifr have already 



been stated in Eqs. (A.5 A.6 1. We are thus left to calculate the 
slightly cumbersome derivatives of multiple partial traces of p 
with respect to the matrix elements pa/j- 



Similarly to the calculation in Appendix A. 1 we will now 
rewrite indices v e {1, . . . , 2^^) of the full Hilbert space in the 
binary representation v = (vi - 1)2'^"^ + (v2 - 1)2'^"^ -H . . . -(- 
2vn^i+v ff = [vi,V2,- ■ ■ ,vai], where V,- € {1,2) for all/ € 1, . . .A^ 
(we will also employ this represention for the indices a,/3, and 
n below). Then, the matrix elements of pt can be written as 

(Pk)iJ — Zj ■ ■ ■ 2j 2j ■ ■ ■ ZjP[Vi,...,Vt_i,l,Vj:+i,...VA.] [Vi,...,Vi_iJ>i+i,...VA.]5 

where i,je {1,2). Hence we have 



d(pk) 



Ji 



dpa/i 



- ^a-,/3, 



3»2j82 



■■5, 



at-lfik- 



OkJ ■ ^/itJ ■ '^ffA-+lA+ 



(A. 13) 
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Inserting this into Eq. ( |A.12| l and working out all Kronecker 6 
symbols, we arrive at 



dy 



d Re i^„ 



N 2 



|iA> 



""X^ZjZj'^K'^I"-"-- 



.«i:-lj,ni-+l. ••■,«« 



■(P^kj)' 



k=\ j=\ 



dy 



dlmij/n 



m 



8 ^ ^ 



k=l j=\ 



(A. 14) 

,«2,...,«i:-i,,/,«i-+i ,...,««] ■ \Pk)ni,j) ■ 

(A.15) 



Appendix B. Derivatives of the function /i([/) 

We will carry out the calculation explicitly only for the deriva- 
tives with respect to the real part of U, but everything works 
analogously for the imaginary part. For the sake of readability, 
we will drop the usual 'ket' notation and write quantum state 
vectors as (/»,■ = \if/j). We write the kth element of if/j as if/ . 

Differentiating Eq. ( [T5| ) with respect to the real part of the 
^ X r matrix U (with matrix elements Uafi) yields 



dh(U) 
dRsUaB 



z 



dpiiU) 



d Re Uat 



miiPiiU)) 



+Pi(U)J] 



dm(i//) 



j=i 



dijj(i^ 



dt 



U) 



^=^, 



dRtU. 



a/J 



where d is the dimension of the total Hilbert space. Note 
that we have specifically emphasized the {/-dependence of the 
Pi and ij/j via Eqs. (J9] [TO] [TTl [T2| . The first derivative in this 
expression is given by 



dpi 



dRtUal: 






dif/'-/^ ' 

dReUap^ 



j ^ 



dReUa/i^ 
2j]R,{,f,y^*6i,^,x'f) 

( I- 



(B.2) 



;■ 



= 25,Q./l^ Re 



Z^-vZ^)'>)^ 



U) 



\l=\ 



= 25,„^^Re({/,yj), 



where we have used in the last step the orthonormality of the xi 
and the fact that Re c* = Re c for any complex number c. 



As for the derivatives of the state vector, we obtain 

dlf/f 1 dipf 1 _3/2 dpi 

— p. 

dUap yfpidUafi 2 ' dUa/l 



-l=Si„ y[Apxf - pf^5iaAp Re [Utp] ^p 



^Su 



.,U) 



l-X)^'-^Re(Ui^)il/\ 
' Pi '^ Pi 



,U) 



(B.3) 



We can now insert Eqs. ( |B.2| l and ( |B.3[ ) into ( |B.l| l. The final 
result (including the derivatives with respect to the imaginary 
part of U from an analogous calculation) reads 

dh{U) 
a Re UaB 



'a 13 

d 

+ PamYj 



J J) dmjip) 






jj) dmjtfr) 



(B.4) 



dh(U) 
— — -— = 2Ai}'im{Uai3)m(il/JU)) 
o im Uap 

(i) dm{^) 






(j) dm(ip) 



■ ' where 




-^.XP - —777 Re(f/./,)^„(f/) 



^Pa{U) Pa{U) 



(B.5) 



(B.6) 
(B.7) 
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